Dual communities in spatial networks

Both human-made and natural supply systems, such as power grids and leaf venation networks, are built to operate reliably under changing external conditions. Many of these spatial networks exhibit community structures. Here, we show that a relatively strong connectivity between the parts of a network can be used to define a different class of communities: dual communities. We demonstrate that traditional and dual communities emerge naturally as two different phases of optimized network structures that are shaped by fluctuations and that they both suppress failure spreading, which underlines their importance in understanding the shape of real-world supply networks.

Supplementary Figure 1. Primal and dual communities emerge naturally in optimal supply networks. a-c, We consider a triangular lattice with two fluctuating sources located at the leftmost and the rightmost corner. The remaining nodes correspond to sinks, whose demand fluctuates following an iid normal distribution. The sources balance the demand in total, but the spatial distribution varies according to a Dirichlet distribution introducing additional fluctuations with variance σ 2 D . We then determine the optimal network structure and edge weights to minimise the average dissipation [1]. We observe a transition from (a) a primal community structure at low fluctuations to (c) a dual community structure at high fluctuations. d-f, Communities can be identified according to the Fiedler vector. The figure shows (d,e) the primal Fiedler vector and (f ) the dual Fiedler vector in a colour map for the networks depicted in panels (a-c), respectively. g, The transition from a primal to a dual community structure can be quantified either by the algebraic connectivity (cf. main text) or by the modularity [2]. The figure shows the modularity Q of the primal community decomposition and the modularity Q * of the dual community decomposition as a function of the variance σ 2 D of the source strength. Supplementary Figure 2. Primal and dual communities emerge naturally in optimal supply networks with three sources. a-c, We consider a triangular lattice with three fluctuating sources located at the leftmost, the upper right and the lower right corner. The remaining nodes correspond to sinks, whose demand fluctuates following an iid normal distribution. The sources balance the demand in total, but the spatial distribution varies according to a Dirichlet distribution introducing additional fluctuations with variance σ 2 D . We then determine the optimal network structure and edge weights to minimise the average dissipation [1]. We observe a transition from (a) a primal community structure at low fluctuations to (c) a dual community structure at high fluctuations. d-f, Communities can be identified according to the Fiedler vector. The figure shows (d,e) the primal Fiedler vector and (f ) the dual Fiedler vector in a colour map for the networks depicted in panels (a-c), respectively. h, The scaling of the corresponding primal (λ2, circles) and dual (λ * 2 , crosses) Fiedler value confirms the transition. With decreasing levels of CO2, generation shifts more and more towards fluctuating renewable power sources with pronounced regional differences: Wind power in predominantly generated along the coast, especially around the North Sea. Solar power is generated predominantly in Southern Europe. Supplementary Figure 5. Definition of a cut-path. a,b, A graph (solid lines) and its weak dual (dashed lines). The edge set of a cut-path S is indicated in red colour. c,d, The cut-path induces a cut of the dual graph into two connected vertex sets V * 1 (blue) and V * 2 (green). The set S * is a minimal cut set in the dual graph. e,f, The faces corresponding to the vertex sets V * 1 and V * 2 form two connected polygons A1 and A2 in the plane.  Figure 6. Different clustering methods in primal and dual graphs Distinct primal and dual communities can be found using both spectral clustering and the Louvain Detection Algorithm. The original (i.e. primal) leaf venation network of Acer platanoides is shown in black in each panel. The resulting hierarchical cuts for three cut levels employing spectral clustering using the primal graph Laplacian and the dual graph Laplacian can be seen in panel a and c. The community structure revealed by the Louvain Detection Algorithm, choosing a resolution parameter ξ = 0.01, using the primal graph and dual graph are shown in panel b and d, respectively.

Supplementary Note 1: Graph duality
In this note, we review the fundamental definitions and properties of dual graphs extending the main text.

Basic definitions
Let G = (V, E) be a plane, undirected graph with N = |V | vertices and L = |E| edges. We label all vertices consecutively as (1, 2, 3, . . . , N ). Edges are either referred to in terms their respective endpoints or by consecutive labels as (1, 2, 3, . . . , L). We begin with unweighted graphs and include edge weights later. For each edge, we fix an orientation -for instance to keep track about the direction of flows as discussed below. The topology of the graph and the orientations are encoded in the node-edge incidence matrix I ∈ R L×N with entries [3] I ℓ,n =    1 if edge ℓ starts at node n, −1 if edge ℓ ends at node n, 0 otherwise.
The Graph Laplacian is defined as This definition coincides with the definition in terms of the adjacency matrix A ∈ R N ×N with entries This property is easy to understand for flow networks: Here, the vector ⃗ v describes a cycle flow which has neither sink or source such that I⃗ v = ⃗ 0. The cycle flows constitute a vector space that coincides with the kernel of the edge incidence matrix I. In plane graphs, the faces provide a distinguished basis of the cycle space. This basis is encoded in the cycle incidence matrix C with entries For convenience, we fix the orientation of each face in the counter-clockwise direction.
We now proceed to the definition of the dual graph G * = (V * , E * ). Two common definitions exist, differing in the treatment of the exterior of the graph.
1. The vertex set V * consists of all faces of the primal graph, including the outer face. Then, the dual graph has |V * | = L − N + 2 vertices.
2. The vertex set V * consists of all faces of the primal graph, excluding the outer face. Then, the dual graph has |V * | = L − N + 1 vertices. This graph is typically referred to as a weak dual.
Two faces (=dual vertices) are adjacent if they share at least one edge. We note that two faces may share more than one edge, which can be included in two different ways. First, we keep every edge separately making the dual graph a multigraph. Alternatively, we can condense the multiedges into ordinary ones, which requires the introduction of edge weights. We will follow the latter definition as our focus is on flow networks which are weighted anyway.
We may now define a Laplacian for the dual graph. Starting with the first definition that includes the outer face, we define, in analogy to Eq. (1), the dual Laplacian as As for the primal graph, we can decompose the Laplacian in terms of the adjacency and the degree matrix, The dual adjacency matrix A * ∈ R (L−N +2)×(L−N +2) has the entries A * f,g = number of edges shared by faces f and g . and the degree matrix In the second definition, we exclude the exterior face. We may again proceed to define the Graph Laplacian via the cycle incidence matrix and write However, this matrix is not a full Laplacian matrix as its row-sums are not guaranteed to vanish. Instead, L * g is a grounded Laplacian, that is obtained from a full Laplacian matrix by removing one row and the corresponding column [4]. Alternatively, we can define the Laplacian via the adjacency matrix as The adjacency matrix A * is a (L − N + 1) × (L − N + 1) matrix has the entries A * f,g = number of edges shared by faces f and g .
as before and the degree matrix Hence, the row-sums of L * i vanish as required for a Laplacian matrix. We note that the two matrices are related as where the diagonal matrix D • summarizes the connectivity of the interior and exterior D • f,f = number of edges shared between f and the exterior face.
Finally, we want to briefly comment on the two alternative definitions of the dual graph and their respective advantages and disadvantages: • The first definition including the exterior face has two closely related advantages making this definition convenient in fundamental analysis: (i) The definitions (3) and (4) of the graph Laplacian coincide. (ii) Every edge is part of exactly two faces. This property forms the basis of MacLane's planarity criterion.
• However, the first definition is rather inconvenient in the analysis of network structures. The cycle basis is overcomplete. The additional exterior face introduces "artificial" non-local connections in the dual, which are especially unsuitable for the definition and analysis of dual communities. For community detection, we will thus work only with the adjacency matrix (7) and the Laplacian (6).
• In the analysis of flow networks we can use both definitions alike. In the first case, we have one degree of freedom that remains to be fixed. A common convention is to "ground" the exterior face -then the first definition reduces to the second one.

Weighted Graphs I
We now consider a weighted plane graph with edge weights denoted as w e ∈ R >0 for all edges e ∈ E. In the analysis of graph duality and flow networks, we can include these weights in two different ways: (i) directly in the incidence matrices and (ii) via separate diagonal weight matrices. We first treat option (i) and focus on the definition of the dual that includes the exterior face.
Subsequently, we define the weighted node-edge incidence matrix I ∈ R L×N as The Graph Laplacian is defined as before with entries To ensure the defining equation for the cycle space, I ⊤ C = 0, the definition of the cycle edge incidence matrix must be generalized to The dual Laplacian matrix is again given by the equivalent definitions (3) and (4), where the adjacency matrix is now given by A * f,g = edges e shared by faces f and g and the diagonal degree matrix is given by D * f,f = g A f,g as usual. We stress that the dual weights are inverse to the primal weights as expressed via Eq. (11).

Weighted Graphs II
A second option is to include the weights via an additional diagonal weight matrix while keeping the unweighted incidence matrixĨ ∈ R L×N with entries 1 if edge ℓ starts at node n, −1 if edge ℓ ends at node n, 0 otherwise.
The Laplacian matrix of the primal graph is then defined by the definition which coincides with Eq. (10) as desired.
We now proceed to the dual, focusing on the second definition that excludes the exterior face. The cycle edge incidence matrixC ∈ R L×(L−N +1) now remains unweighted with components The dual grounded Laplacian (5) is recovered via the definition Comparing to the definition of the primal Laplacian (14), we again see that dual weights are inverse to primal weights.

Supplementary Note 2: Algebraic and topological connectivity
The correspondence of algebraic and topological connectivity is well established for primal graphs. Here, we briefly review this correspondence and extend it to the dual graph.

Primal Graph
Let G = (V, E, W) be a weighted connected graph and L its Laplacian. The Laplacian has one vanishing eigenvalue λ 1 = 0 with corresponding eigenvector ⃗ v 1 ∝ ⃗ 1 = (1, 1, . . . , 1) ⊤ . Since L is positive semi-definite [3], we can write the eigenvalues in an ordered way The second eigenvalue λ 2 , the Fiedler value, provides a measure of the algebraic connectivity of a graph. If it is small, we have a pronounced community structure. We now relate the algebraic connectivity λ 2 to the topological connectivity across a cut-set of the graph. To this end, we first express λ 2 via the variational principle as We can then find an upper bound by choosing an appropriate trial vector for ⃗ v. To this end, we decompose the vertex set into two parts V 1 and V 2 = V \V 1 and define a trial vector with entries where N 1 = |V 1 | and N 2 = |V 2 |. Then we have We note that the expression (20) is not just an upper bound but gives λ 2 up to linear order in the topological connectivity. To make this more precise, assume that we start from a graph that is disconnected into two parts V 1 and V 2 = V \V 2 . We label the vertices as V 1 = {1, . . . , N 1 } and V 2 = {N 1 + 1, . . . , N 2 + N 2 }. The adjacency matrix then assumes a block form 2 ∈ R N2×N2 . In this limit we have λ (0) 2 = 0. Now add some weak connections between the two modules of the graph. Then, the adjacency matrix can be written as We now treat A (1) as a small perturbation and apply Rayleigh-Schrödinger perturbation theory. To zeroth order we have λ coincides with the trial vector defined by Eq. (19). To linear order, the Fiedler value is then given by

Dual Graph
We can now generalize the above treatment to the dual graph. As discussed above, we base the analysis of communities and connectivity on the dual graph that excludes the exterior face. The Laplacian of the dual graph is defined by (6) or its weighted counterpart. In particular, we have where the entries of the adjancency matrix are given by A * f,g = edges e shared by faces f and g and D * f f = g A * f g . As before, we want to define a decomposition of the graph and use the connectivity at the boundary to derive a bound for the algebraic connectivity. However, we now need the weights along the boundary, not across the boundary.
To make this precise, we first define a cut-path that we will use instead of a cut-set. Definition 1. Let G = (V, E) be a connected plane graph. The outer boundary B is defined as the cyclic path consisting of the edges adjacent to the exterior face and the connecting vertices.
that includes no edges from the outer boundary B and that satisfies 1. either the path is cyclic and contains at most one vertex from the boundary B 2. or the path is acyclic and its start and end are on the boundary, v 1 , v n+1 ∈ B.
We will now show that such a cut-path p indeed cuts the dual graph into two parts. Then we show that the sum of weights w −1 e along the cut-path provides an upper bound for the algebraic connectivity of Fiedler value of the dual. Lemma 1. Let G = (V, E) be connected plane graph and G * = (V * , E * ) its dual graph. We consider a cut of the dual G * , i.e. a decomposition of the vertex set V * such that Then the following statements are equivalent 1. The two dual vertex sets V * 1 and V * 2 are connected in G *

The cut set
is minimal.
3. The dual of the cut set S = {e ∈ E|e is dual to an e * ∈ S * } = {e ∈ E|e is adjacent to a face f ∈ V * 1 and a face g ∈ V * 2 }.
is the edge set of a cut-path.
Proof. The equivalence of (1) and (2), i.e. the correspondence of connectedness and minimality of the cut set, is a standard result in graph theory (see, e.g., [5,Lemma 1.9.4.]). The equivalence of (1,2) and (3) follows from the geometry of the plane embedding: A vertex v ∈ V corresponds to a point in the plane R 2 , an edge to an arc and a face to a closed polygon. An exact proof for dual graphs that include the exterior face can be found in many text books of graph theory, see e.g. [5, Proposition 4.6.1.]. Here we briefly sketch how to extend this result to dual graphs that exclude the outer face.
(1) ⇒ (3) Consider the embedding of the graph in the plane R 2 . Then all faces of the graph define a connected area A ⊂ R 2 , whose boundary ∂A corresponds to the edges in the boundary set B. Similar, the faces in V * 1,2 correspond to areas A 1,2 ⊂ R such that A 1 ∪ A 2 = A. As V * 1,2 are connected by assumption, so are the areas A 1,2 . Examples of the possible geometries are shown in Supplementary Figure 5.
The boundary between the two areas ∂A 1,2 = A 1 ∩ A 2 is a connected polygonal chain. It contains at most two points in the boundary ∂A of the area A, at most one if it is closed, cf. Supplementary Figure 5. Now the boundary coincides with the set S = {e ∈ E|e is adjacent to a face f ∈ V * 1 and a face g ∈ V * 2 }.
which thus inherits the properties of ∂A 1,2 : (i) As ∂A 1,2 is a connected polygonal chain, S is the edge set of a path.
(ii) As ∂A 1,2 contains at most two points in ∂A, S contains at most two vertices from B, at most one if the path is cyclic.
(3) ⇒ (1) Let S be the edge set of a cut-path p. We denote this cut-path as p = (v 1 , e 1 , v 2 , . . . , v n+1 ) and distinguish two cases. (a) If the cut-path p is acyclic, its start v 1 and end v n+1 are distinct and both lie on the boundary B. We can thus write the boundary as and define two cyclic pathsp In the plane embedding,p 1 andp 2 correspond to two closed polygon chains enclosing two polygons A 1,2 . As thep 1,2 are cyclic paths, the polygons A 1,2 are connected. We can again identify the areas A 1,2 with the dual vertex sets V * 1,2 and conclude the following statements. (i) Asp 1,2 are cyclic paths, the polygons A 1,2 are simply connected areas. Hence, the dual vertex sets V * 1,2 are connected. (ii) Asp 1,2 contain the boundary B, we have A = A 1 ∪ A 2 and hence V * = V * 1 ∪ V * 2 . (iii) As the pathsp 1,2 do not intersect, the two polygons are disjoint, A 1 ∩ A 2 = ∅, and so are the dual vertex sets, V * 1 ∩ V * 2 = ∅. (b) If the cut-path p is cyclic, its start and end coincide and lie on the boundary v 1 = v n+1 ∈ B. In the plane embedding, the cut-path p corresponds to a closed polygon chain enclosing a polygon A 1 . As p is a cyclic path, A 1 is connected. Notably, A 1 touches the boundary in exactly one point corresponding to the node v 1 = v n+1 according to the definition of a cut-path. The area A 2 = A\A 1 must also be connected, otherwise A 1 would touch the boundary in more than one point. We can now identify the areas A 1,2 with the dual vertex sets V * 1,2 and conclude that both V * 1 and V * 2 are connected and satisfy V * Corollary 1. Let G = (V, E, W ) be a weighted plane graph and be a cut-path. Then the algebraic connectivity of the dual graph G * satisfies where E(p) = {e 1 , . . . , e n } is the set of edges in the cut-path and N * 1,2 = |V * 1,2 | Proof. In full analogy to the analysis of the primal graph (cf. (20)), we find the inequality where V * 1 and V * 2 = V * \V * 2 is a decomposition of the vertex set of the dual and N * 1,2 = |V * 1,2 |. The entries of the dual adjacency matrix are given by Eq. (21) such that we have edges e shared by faces c and d where the last equality follows from Lemma 1.

Supplementary Note 3: Linear flow networks
Linear flow networks describe a generic model for many types of supply networks including AC power grids [6][7][8], DC electric circuits [9][10][11], hydraulic networks [12,13], and vascular networks of plants [14]. Assume that the underlying network is given by a graph G = (E, V ) with N = |V | nodes and M = |E| edges. Then we assign a potential θ n ∈ R to each node n ∈ V . This potential has the following interpretation: In AC power grids in the linear power flow approximation, it denotes the voltage phase angle whereas in DC electric circuits, it refers to the voltage level at node n. Finally, in hydraulic or vascular networks, θ n denotes the pressure at node n. Then the flow F ℓ along a link ℓ = (m, n) is given by where w ℓ is the weight of the link ℓ. In AC power grids in the linear power flow approximation, w ℓ is proportional to the link's susceptance, in resistor networks it is given by the link's conductance and in a hydraulic or vascular network it depends on the geometry of the pipe or vein. Note that all these networks are undirected, i.e. flow is possible in both directions alike, such that w m,n = w n,m . Since the flow (28) depends only on the potential drop, the potentials θ n are only determined up to a constant phase shift applied to all nodes. This problem may be solved by assigning a reference potential to a preselected slack node n, e.g. setting θ n := 0. Now assume that there is an inflow P m at every node m ∈ V . The value of P m is positive if a current, power, or fluid is injected to the node and negative if it is withdrawn from the node. Here and in the following, we assume that the in-and outflows sum to zero, N i P i = 0, and call them 'balanced'. The flows F ℓ along the edges ℓ ∈ E may then be determined by using the continuity equation, also referred to as Kirchhoff's current law (KCL) ℓ∈EĨ ℓ,n F ℓ = P n , ∀ n ∈ V.
Here,Ĩ ℓ,n are the entries of the graph's unweighted edge-node incidence matrix (13). The potentials θ n which fulfil the continuity equation (29) Defining a vector of flows ⃗ F = (F 1 , .., F M ) ⊤ ∈ R M , a vector of potentials ⃗ θ = (θ 1 , ..., θ N ) ⊤ ∈ R N and a vector of inflows ⃗ P = (P 1 , ..., P N ) ⊤ ∈ R N , we can write these relationships in a more compact form. Kirchhoff's current law (29) in vectorized form reads asĨ In addition to that, we can write the relationship between flows and potentials in vectorized form, Here, W = diag(w 1 , ..., w L ) ∈ R M ×M is a diagonal matrix summarizing the edge weights, cf. Eq. 12. Now we can put together the last two equations to arrive at a discrete Poisson equation for the nodal potentials Here, L =Ĩ ⊤ WĨ is the Laplacian matrix.
Quantifying the impact of perturbations -primal graph We now use the notation developed in the last section to study the effect of perturbations, in particular links failures, and derive the expression for the sensitivity factor η i,k,ℓ given in the main text. Assume that the inflow at one node i is increased by an amount ∆P and decreased at another node k by the same amount. The inflow vector changes as ⃗ P → ⃗ P ′ = ⃗ P + ∆ ⃗ θ with where ⃗ e i denotes the ith standard unit vector. As a consequence, the nodal potentials changes as ⃗ θ → ⃗ θ ′ = ⃗ θ + ∆ ⃗ θ. Both the new and the old potentials have to fulfil the Poisson equation (32), Subtracting the two equations, we arrive at an explicit equation for the change in nodal potentials Note that the Laplacian matrix of a connected graph always has one zero eigenvalue [3]. Thus, we cannot simply invert this equation to calculate the change in the nodal potentials ∆ ⃗ θ. This problem is typically solved by making use of the Moore-Penrose pseudoinverse of the Laplacian matrix L † which has properties similar to the matrix inverse [15]. Now we can use Kirchhoff's current law (29) to calculate the resulting changes in the flows Finally, we can read off the change of the flow ∆F ℓ on a line ℓ. We define an indicator vector ⃗ l ℓ as the ℓth standard unit vector in ⃗ R M and obtain where we replaced the unweighted incidence matrixĨ by the weighted matrix I. Then the sensitivity factor reads Note that the sensitivity factor is purely topological, i.e. it may be calculated only based on the network topology and weights and does not depend on the in-and outflows. If the two nodes i and k are the terminal nodes of an edge p = (i, k), we can identify and the sensitivity factor reads Quantifying the impact of perturbations -dual graph The discussion of perturbations in the in-and outflows developed in the last section was recently extended to the dual graph [16,17]. We will briefly cover this here and derive the dual expression for the sensitivity factor η. As before, we assume that the inflow at one node i is increased by an amount ∆P and decreased at another node k by the same amount such ∆ ⃗ P = ∆P (⃗ e i − ⃗ e k ). By subtracting Kirchhoff's current law (29) for the situation before and after the change in inflow, we arrive atĨ which leaves us with an underdetermined system of equations: We have M unknown flow changes ∆F ℓ , but only N − 1 equations as one is redundant. The general solution can thus be written as where ∆ ⃗ F part is one particular solution of the equation and ∆ ⃗ F hom is an arbitrary solution of the homogeneous equation, i.e. an arbitrary solution ofĨ ⊤ ∆ ⃗ F hom = ⃗ 0. We can now use graph duality to parametrize the homogeneous solutions and to single out the physically correct solution.
The kernel of the node-edge incidence matrix is spanned by the cycle-edge incidence matrix where we here use the unweighted matrix (15). Hence, we can write the homogeneous solutions as is a vector of cycle flows -one for each cycle. The physically correct values of the cycle flows are then determined by Kirchhoff's voltage law (Eq. 30) which may be written compactly as where we used that ∆θ n − ∆θ m = ∆F nm /w nm . Inserting equations (35) and (36), we then obtaiñ with the dual Laplacian L * g =C ⊤ W −1C . We thus found a discrete Poisson equation for the cycle flows ⃗ f with the grounded dual Laplacian L * g in direct correspondence to its primal counterpart (32). We can now proceed as for the primal graph: Inverting the discrete Poisson equation for the cycle flows and plugging the result into Eq. (35), we thus arrive at Here, 1 M denotes the identity matrix of dimensions M × M . As a last step, we need to determine a particular solution ∆ ⃗ F part . This can be accomplished as follows [17]: Let T ik ∈ R M be a vector determining an (arbitrary) path between the node i with inflow ∆P and node k with the same outflow. The entries of T ik are given as follows Then a particular solution is given by ∆ ⃗ F part = T ik ∆P , and we can plug this into the above equation If nodes i and k are the two terminal ends of an edge e = (i, k) ̸ = ℓ, the path vector may be identified with the edge's indicator vector, T ik = ⃗ l e . Then we can simplify the above expression for any edge ℓ ̸ = e to Notably, the first expression in Eq. (38) vanishes since ⃗ l ⊤ ℓ ⃗ l e = δ ℓe . Thus, we can identify the sensitivity factor η i,k,ℓ as

Modelling link failures
The sensitivity factor η i,k,ℓ may also be used to describe the failure of links. Assume that a link e fails, thus loosing its ability to carry flow and setting the weight to zero, w e = 0. Instead of removing the link from the network, we can find an analogous description based on the sensitivity factor [6,18].
Assume that the link e = (e 1 , e 2 ) carries a flow F e before the outage. Now we model the removal of the link by an in-and outflow at the two ends of the link: To this end, we assume that the line was disconnected from the network and consider a fictitious flowF e which is the result of a (fictitious) inflow ∆P =F e at the starting node e 1 and an outflow of the same amount at the terminal node e 2 of link e. On the other hand, we can also calculate the flowF e flowing on line e after the injection by a self-consistency argument: It can be calculated using the sensitivity factor F e = F e + η e1,e2,e ∆P = F e + η e1,e2,eFe ⇒ ∆P =F e = F e (1 − η e1,e2,e ) −1 .
Now we can calculate the change in the flow on another link ℓ due to the failure of link e as ∆F ℓ = η e1,e2,ℓ ∆P = η e1,e2,ℓ 1 − η e1,e2,e F e .
We can thus calculate the flow change ∆F ℓ on any link ℓ due to the failure of another link e based on the initial flow on link e and the sensitivity factor η. The expression η e1,e2,ℓ (1 − η e1,e2,e ) −1 is known as Line Outage Distribution Factor in power system security analysis [6]. Thus, both changes in the inflows and link failures can be captured by the sensitivity factor η.
To quantify the effect of (dual) communities on network robustness in linear flow networks we defined the ratio of flow changes in the main text (cf. Ref. [19]), Importantly, the ratio can be used to quantify both, the effect of communities on changes in the inflow patterns and on failure spreading. For an inflow and outflow of ∆P at the two terminal ends of the link e, e 1 and e 2 , respectively, the ratio is calculated as On the other hand, if we instead calculate the ratio for the failure of a link e with initial flow F e , we arrive at Thus, in both cases, the ratio is determined by the sensitivity factor η which is in turn governed by the Pseudo-inverse of the Laplacian matrix L † or the Pseudo-inverse of the dual Laplacian (L * g ) † when formulating the problem in the dual graph.

Connectivity structure determines network response to link failures
In this section, we will demonstrate why a weak connection between two components of a network limits the flow changes in one component when a link in the other one fails. Our analysis follows the approach towards perturbation spreading used in Manik et al. [20,21] that is based on Rayleigh-Schrödinger perturbation theory [22].
Consider a connected graph G = (E, V ) with N nodes and M edges consisting of two subgraphs G 1 = (E 1 , V 1 ) and G 2 = (E 2 , V 2 ) with n 1 nodes and n 2 = N − n 1 nodes, respectively, that are mutually weakly connected. Here, a weak connection between the two subgraphs may either be realized through a (relatively) weak number of links in case of an unweighted graph or the overall weight of the connections between the two subgraphs being small [5]. We then sort the graph's vertices V in such a way that the first n 1 vertices belong to the first subgraph G 1 and the other n 2 vertices belong to the second one G 2 . We may regard the weak connections between the to subgraphs as a perturbation to the graph in which the two subgraphs are disconnected. In terms of the graph Laplacian L, we thus write Here, L 0 is the graph Laplacian for the graph when disconnecting the two subgraphs, expressed in terms of their Laplacian matrices L 1 ∈ R n1×n1 and L 2 ∈ R n2×n2 andL is the perturbation matrix with the diagonal degree matrices D 12 and D 21 denoting the degree for the graph connecting the two subgraphs. A 12 is the adjacency matrix for this graph that is assumed to be relatively sparse, thus indicating the weak inter-subgraph connections.
To examine the effect of link failures, we need to study the pseudoinverse of this matrix. Therefore, we denote by X = L † the Moore-Penrose pseudoinverse of the overall graph Laplacian and by X 1 = L † 1 and X 2 = L † 2 the Moore-Penrose pseudoinverses of the subgraphs' Laplacian matrices L 1 and L 2 , respectively. For the inverse X 0 of the unperturbed Laplacian L 0 , we then get In order to calculate the matrix inverse X, we can expand this matrix using the Neumann series (see e.g. Ref. [23]). We note that some issues may arise in this perturbative treatment and the series expansion, as the Moore-Penrose pseudo-inverse is not continuous in general. However, these this does not apply here, since we can easily circumvent the formal problems by fixing the phase of an arbitrarily chosen slack node n as θ n ≡ 0 and removing the n-th row from the linear set of equations. The resulting set of equations is of full rank and admits a unique solution. The corresponding matrices are again grounded Laplacians [4] which are invertible. The matrix inverse is continuous and admits a series expansion as in Equation (39). For the sake of notational simplicity and coherence, we will however stick to the ordinary Laplacians. The matrix inverse X then reads where we inserted the Neumann series in the last step. We can thus approximate the matrix inverse of the graph Laplacian as where O(L 2 ) denotes terms of at least order two in the perturbation matrixL. Now we can use these expressions to calculate a first order approximation for the sensitivity factor η i,k,ℓ in the weakly connected limit. Assume we are monitoring the flow changes on line ℓ with indicator vector ⃗ l ℓ and define ⃗ ν ℓ = I · ⃗ l ℓ = ⃗ e ℓ1 − ⃗ e ℓ2 as a result of a power transfer along line k = (k 1 , k 2 ) with indicator vector ⃗ l k and ⃗ ν k = I · ⃗ l k . We can then calculate the sensitivity factor as η k1,k2,ℓ = w ℓ ⃗ ν ⊤ ℓ X⃗ ν k . Now we distinguish two cases. First assume that ℓ and k are contained in the same subgraph, say G 1 . In this case, we can write with slight abuse of notation ⃗ ν ℓ = (⃗ ν ℓ , ⃗ 0 n2 ) ⊤ and ⃗ ν k = (⃗ ν k , ⃗ 0 n2 ) ⊤ , where⃗ ν l ∈ R n1 denotes the projection of the vector onto the subspace of vertices the first subgraph G 1 and the whole vector ⃗ ν ℓ is still understood as a vector in R N . In this case, we may write the sensitivity factor as Now consider the case where ℓ and k are contained in different modules of the network. In this case, we may write ⃗ ν ℓ = (⃗ ν ℓ , ⃗ 0 n2 ) ⊤ and ⃗ ν k = ( ⃗ 0 n1 ,⃗ ν k ) ⊤ and calculate the sensitivity factor as We thus observe that the leading order contribution of the perturbation matrix to the sensitivity factor isL 0 if both links are contained in the same module andL 1 if they are contained in different modules. Now we compare this to the scaling of Laplacian eigenvalues with the perturbation matrix. The first eigenvector of L 0 is given by the constant shift ⃗ v 1 = N −1/2 ⃗ 1 N and has eigenalue λ 1 = 0. In the weakly connected limit considering only L 0 , the Fiedler eigenvalue vanishes as well λ (0) 2 = 0, and has an associated eigenvector [20] ⃗ v Here, we use the superscript (0) to denote the eigenvector and eigenvalue in the unperturbed case. A first order estimate for the Fiedler value may thus be calculated by using Rayleigh Schrödinger perturbation theory as For weakly connected graphs, we thus expect the sensitivity factor η ℓ1,ℓ2,k to scale with the connectivity of the graph in the same way as the Fiedler value of the graph if ℓ and k lie in different communities due to the fact that λ 2 ∝L and expect it to be to leading order independent of the Fiedler value if ℓ and k are in the same community. We start by a review of the original single-sink model as formulated by Corson [1]. Consider a linear flow network on a graph G with N nodes and M edges summarized in the node set V and edge set E. Then we choose one vertex to be the source of the network while all the others are sinks. We assume that the outflow at the sinks fluctuates in time, which is modelled by uncorrelated, iid Gaussian variables P i ∼ N (µ, σ), i ∈ {2, ..., N }. Sources and sinks have to balance at any point in time such that This equation immediately yields the statistical properties of the source as a consequence of the statistics of the sinks. We label the nodes such that the source has the index i = 1 and obtain where ⟨·⟩ denotes the mean over different realizations of the Gaussian variables. The network structure described by the edge weights w ℓ is then optimized to minimize the average dissipation Here, ⟨F 2 ℓ ⟩ is the second moment of the flows which are determined by the statistics of the sources and sinks by virtue of Kirchhoff's equations. However, we assume that the budget for constructing or strengthening edges is limited leading to a resource constraint ℓ∈E w γ ℓ ≤ 1.
Here, γ is a cost parameter that controls how expensive it is to increase the weight of an edge. The main challenge when performing the minimization is the interdependence between link weights and flows which cannot be varied independently. Solving the constrained optimization problem with using the method of Lagrange multipliers yields two conditions: First, the flows are determined by the edge weights w ℓ and the inflows P i via Kirchhoff's equations. Second, the optimal weights are related to the flows as Corson [1] and Katifori et al. [14] independently came up with a procedure to solve this problem in an iterative self-consistent manor starting from a random initial guess for the edge weights. Then one out the following steps repeatedly until the procedure converges: 1. The moments of the flows are determined via Kirchhoff's equations.
Notably, the system can admit several local minima, and it depends on the initial guess, to which minimum the iteration converges. Now we extend this set-up to a network with multiple sources. Assume that we again label the nodes such that the first N s nodes are sources and the remaining N − N s nodes are sinks, which are still uncorrelated, iid Gaussian variables. However, when considering more than one source, N s > 1, the distribution of the sources is not completely determined by Equation (41), but has additional degrees of freedom. We use this degree of freedom to put additive Dirichlet noise X i ∼ Dir(α) on the sources. In particular, we assume that where K ∈ R is a scaling parameter. Importantly, the Dirichlet variables sum to unity for all realizations, The parameter α determines the shape of the dirichlet distribution and tunes the variability of the sources. For α = 1, the distribution is flat such that the individual X i fluctuate strongly. For α ≫ 1, the distribution strongly peaks around the mean such that fluctuations of the individual X i are small. The moments of the Dirichlet random variables are given by s (N s α + 1) Note that this results in a vanishing sum over the second moments by virtue of the definition of the Dirichlet parameter and the first equation. Now we can shift to new random variables Y i = K 1 Ns − X i with zero mean, ⟨Y i ⟩ = 0, and the following second moments such that the sum over the second moment still vanishes Ns i=1 ⟨Y i Y j ⟩ = 0. With these results we can compute the second moments of the inflows ⟨P i P j ⟩. Given the network topology and the edge weights w ℓ , we can further compute the moments of the edge flows ⟨F 2 ℓ ⟩ using Kirchhoff's equations. Then one can proceed as in the original model as outlined above. Starting from an initial guess of the edge weights, we iterate the following two steps until convergence: 1. The moments of the flows ⟨F 2 ℓ ⟩ are determined from the moments ⟨P i P j ⟩ via Kirchhoff's equations.

Supplementary Note 5: Clustering in primal and dual graphs
While we used spectral clustering to find communities in both primal and dual graphs to enable an analytical treatment in the main manuscript, there are various different algorithms to determine the community structure in a given network. A commonly used and efficient algorithm is the Louvain Community Detection Algorithm, which is based on finding a community partition by iteratively optimizing the modularity of this partition. We are using the algorithm implemented in Networkx-2.8.6 [27], which uses a combination of [28] and [29].
The algorithm consists of two steps that are run iteratively. Starting from a partition in which every node is in its own community, the change in modularity ∆Q for removing an isolated node i from their lonely community and moving them to a neighboring community C is calculated using where m is the total number of edges in the graph, k i,in is the sum of the weights of the edges connecting node i to the community C, and S tot is the sum of the weights of all edges connected nodes in C. A resolution parameter ξ smaller than 1 favors larger communities, while ξ larger than 1 favors smaller communities. The first step continues moving isolated nodes to neighboring communities until no further increase in modularity gain can be achieved. In the second step, the previously found communities are interpreted as nodes of a new network, and the weights of the new edges that connect the new nodes are given be the sum of the edges that connected the different communities. Subsequently, the first step can be executed again using this new network to find larger communities. These two steps are run iteratively until no increase in modularity that is larger than a small threshold is achieved. Choosing a resolution parameter of ξ = 0.01 and running the clustering algorithm for the primal and dual graph of the Acer platanoides leaf, we find the community structure of the given graphs and compare them with the ones found by hierarchical spectral clustering. The results are shown in Fig. 6. Note, the set of edges that separate the communities are the edges connecting the different communities in the case of primal graphs, while they are along a path in the case of dual communities. In both cases the primal (i.e. original) graphs are shown in the background in each panel in Fig. 6. The separation along veins is even more pronounced and accurate when using the Louvain method, indicating that the description as dual graphs reveals important structural features of the topology of leaf venation networks.